######################################################################################
######################################################################################
#Replication code for models presented in:

#Buntaine, Mark T., Bradley C. Parks, and Benjamin P. Buch
#Aiming At the Wrong Targets: The Domestic Consequences of International Efforts to Build Institutions
#International Studies Quarterly

#Mark Buntaine contact (as of September 2016): buntaine@bren.ucsb.edu / 805.893.4075
#Bradley Parks contact (as of September 2016): bcpark@wm.edu / 757.221.2463
#Benjamin Buch contact (as of September 2016): bpbuch@stanford.edu

#Current for R Version 3.3.1 (version "Bug in Your Hair")
######################################################################################
######################################################################################

##########################################################################################
###Setting Up
##########################################################################################
#install.packages("devtools")
#install.packages("ez")
#install.packages("Hmisc")
#install.packages("ROCR")
#install.packages("lme4")
#install.packages("foreign")
#install.packages("arm")
library(devtools) #Version 1.12.0
library(ez) #Version 4.3
library(Hmisc) #Version 3.17-4
library(ROCR) #Version 1.0-7
library(lme4) #Version 1.1-12
library(foreign) #Version 0.8-66
library(arm) #Version 1.9-1

data <- read.csv("~/Dropbox/PIDs/Round4/Data_InSample.csv", stringsAsFactors=FALSE)
data2 <- read.csv("~/Dropbox/PIDs/Round4/Data_OutofSample.csv", stringsAsFactors=FALSE)

#Re-leveling the environmental ministry data
data$env.min.combined <- as.factor(data$env.min.combined)
data$env.min.combined <- relevel(data$env.min.combined, ref = "None")

data2$env.min.combined <- as.factor(data2$env.min.combined)
data2$env.min.combined <- relevel(data2$env.min.combined, ref = "None")

data$op.cat.I <- ifelse(data$WB.operational.cat=="I",1,0) #Creating indicator for WB operational category I (IDA only)
data2$op.cat.I <- ifelse(data2$WB.operational.cat=="I",1,0) #Creating indicator for WB operational category I (IDA only)

########################################################################
###Table 2
########################################################################

m1 <- glmer(result ~ ida.prop + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
summary(m1)

m2 <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
summary(m2)

m3 <- glmer(result ~ op.cat.I + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
summary(m3)

########################################################################
###Figure 2
########################################################################

sim.gate1 <- sim(m1,n.sims=10000)
fixef1 <- fixef(sim.gate1)
sim.gate2 <- sim(m2,n.sims=10000)
fixef2 <- fixef(sim.gate2)
sim.gate3 <- sim(m3,n.sims=10000)
fixef3 <- fixef(sim.gate3)

#pdf(file="~/Dropbox/PIDs/Round4/TargetChoice_Figure1.pdf", width=7.5,height=7.5)
par(mfrow=c(2,2),oma=c(0,0,0,0.1),mar=c(4.3,4.3,4.3,0.3))

sim1.base <- invlogit(fixef1[,1] + fixef1[,2]*0 + fixef1[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef1[,6]*mean(data$approval_year.count,na.rm=T) + fixef1[,7]*0)
sim1.IV <- invlogit(fixef1[,1] + fixef1[,2]*1 + fixef1[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef1[,6]*mean(data$approval_year.count,na.rm=T) + fixef1[,7]*0)
sim1.diff <- sim1.IV-sim1.base
plot(density(sim1.diff),lwd=3,xlim=c(-0.4,0.1),ylim=c(0,5.5),
     xlab="Change in Probability of Function Target \n(Baseline Probability = 0.56)", cex.lab=1.2,
     main="M1. IBRD Project to IDA Project",cex.main=1.1)
abline(v=0,lty=2,lwd=5,col="red")
text(x=-0.05,y=5,paste("99.3%"),cex=1.2)
text(x=0.05,y=5,paste("0.7%"),cex=1.2)

sim2.base <- invlogit(fixef2[,1] + fixef2[,2]*0 + fixef2[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef2[,6]*mean(data$approval_year.count,na.rm=T) + fixef2[,7]*0)
sim2.IV <- invlogit(fixef2[,1] + fixef2[,2]*1 + fixef2[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef2[,6]*mean(data$approval_year.count,na.rm=T) + fixef2[,7]*0)
sim2.diff <- sim2.IV-sim2.base
plot(density(sim2.diff),lwd=3,xlim=c(-0.4,0.1),ylim=c(0,5.5),
     xlab="Change in Probability of Function Target \n(Baseline Probability = 0.56)", cex.lab=1.2,
     main="M2. IBRD Portfolio to IDA Portfolio",cex.main=1.1)
abline(v=0,lty=2,lwd=5,col="red")
text(x=-0.05,y=5,paste("95.8%"),cex=1.2)
text(x=0.05,y=5,paste("4.2%"),cex=1.2)

sim3.base <- invlogit(fixef3[,1] + fixef3[,2]*0 + fixef3[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef3[,6]*mean(data$approval_year.count,na.rm=T) + fixef3[,7]*0)
sim3.IV <- invlogit(fixef3[,1] + fixef3[,2]*1 + fixef3[,3]*mean(data$NR.total.rents.ay.lag1,na.rm=T) + fixef3[,6]*mean(data$approval_year.count,na.rm=T) + fixef3[,7]*0)
sim3.diff <- sim3.IV-sim3.base
plot(density(sim3.diff),lwd=3,xlim=c(-0.35,0.1),ylim=c(0,6),
     xlab="Change in Probability of Function Target \n(Baseline Probability = 0.54)", cex.lab=1.2,
     main="M3. IBRD Eligible to IDA-Only Category",cex.main=1.1)
abline(v=0,lty=2,lwd=5,col="red")
text(x=-0.05,y=5,paste("98.0%"),cex=1.2)
text(x=0.05,y=5,paste("2.0%"),cex=1.2)
#dev.off()

########################################################################
###Table 3
########################################################################

#Creating subset of data with common support on per capita income
data.band <- subset(data, GDP.k.pc.ay.lag1>=0.3985 & GDP.k.pc.ay.lag1<1.4710)

m4 <- glmer(result ~ ida.prop + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.band, control=glmerControl(optimizer="bobyqa"))
summary(m4)

m5 <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.band, control=glmerControl(optimizer="bobyqa"))
summary(m5)

#Note: rescaling NR.total.rents.ay.lag1 addresses convergence issue with m5, leading to same estimates
data.band$NR.total.rents.ay.lag1.prop <- data.band$NR.total.rents.ay.lag1/100
m5.rescaled <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1.prop + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.band, control=glmerControl(optimizer="bobyqa"))
summary(m5.rescaled)

m6 <- glmer(result ~ op.cat.I + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.band, control=glmerControl(optimizer="bobyqa"))
summary(m6)

#Creating subset of data within margin of error on operational categorization based on historical revisions to the Atlas data
data.catI.error <- subset(data, cutoff.diff>=-375 & cutoff.diff<=175)

m7 <- glmer(result ~ ida.prop + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.catI.error)
summary(m7)

m8 <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.catI.error, control=glmerControl(optimizer="bobyqa"))
summary(m8)

m9 <- glmer(result ~ op.cat.I + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.catI.error)
summary(m9)

#Creating subsets based on categorization, to examine the subgroup effects of income
data.i <- subset(data, op.cat.I==1)
data.not.i <- subset(data, op.cat.I==0)

m10 <- glmer(result ~ GDP.k.pc.ay.lag1 + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.i, control=glmerControl(optimizer="bobyqa"))
summary(m10)

m11 <- glmer(result ~ GDP.k.pc.ay.lag1 + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=data.not.i, control=glmerControl(optimizer="bobyqa"))
summary(m11)
  
########################################################################
###Figure 3
########################################################################

#Figure 3 Model 1
pred <- predict(m1, newdata = data2, type = "response", allow.new.levels = TRUE)
predict <- prediction(pred, data2$result)
perf <- performance(predict,"tpr", 'fpr')

##Figure 3 Model 2
pred2 <- predict(m2, newdata = data2, type = "response", allow.new.levels = TRUE)
predict2 <- prediction(pred2, data2$result)
perf2 <- performance(predict2,"tpr", "fpr")

##Figure 3 Model 3
pred3 <- predict(m3, newdata = data2, type = "response", allow.new.levels = TRUE)
predict3 <- prediction(pred3, data2$result)
perf3 <- performance(predict3,"tpr","fpr")
#performance(predict3,"tpr","fpr", measure="auc") #Example code to get AUC value

#pdf("Round4_Figure2.pdf",width=6.5,height=6.5)
par(mfrow=c(2,2), mar=c(2.5,2.5,2.5,2.5)+0.1, mgp=c(1.5,.5,0), cex.axis = .6)
plot(perf,main="M1. IDA Project Model Predictive Power", lwd=2, cex.lab = .7, cex.main = .8)
abline(0,1, lty =3)
text(0.3,0.2,"Total Area Under Curve: 0.654",adj=c(0,1), cex = .7)
plot(perf2,main="M2. IDA Portfolio Model Predictive Power", lwd=2, cex.lab = .7, cex.main = .8)
abline(0,1, lty = 3)
text(0.3,0.2,"Total Area Under Curve: 0.650",adj=c(0,1), cex = .7)
plot(perf3,main="M3. IDA-Only Category Predictive Power", lwd=2, cex.lab = .7, cex.main = .8)
abline(0,1, lty = 3)
text(0.3,0.2,"Total Area Under Curve: 0.662",adj=c(0,1), cex = .7)
#dev.off()

########################################################################
###Figure 4
########################################################################

#Treeplot function to plot predictive power of individual covariates
treeplot2 <- function(model,DV,IVs,IVnames, newdv, olddat, newdat,x1,y1,x2,downshift=0.002,rightshift=0.015) {
  # first calculate the AUCs by deleting one variable at a time:
  AUCs<-rep(NA,5)
  names(AUCs)<-names(IVs)
  vars <- names(IVs)
  for (i in 1:ncol(IVs)) {
    vars.current <- paste(names(IVs)[-i],collapse=" + ")
    formula.current <- paste("result ~ ", vars.current, " + (1|level) + (1|country)", sep="")
    model.current<-glmer(formula=formula.current, family = binomial, data=olddat, control=glmerControl(optimizer="bobyqa"))
    arg <- plogis(predict(model.current, newdata = newdat, type = "response", allow.new.levels = TRUE))
    AUCs[i]<-somers2(arg, newdv)[1]
  } # close i loop
  
  # now plot the points and lines:
  for (i in 1:length(vars)) {
    segments(x1,y1,x2,AUCs[i])
    points(x2,AUCs[i],cex=1,pch=21,bg="green")
  }
  # now add the text:
  textheight<-AUCs[order(AUCs)]
  orderednames<-IVnames[order(AUCs)]
  for (i in (ncol(IVs)-1):1){
    if (textheight[i]>textheight[i+1]-downshift)
      textheight[i]<-textheight[i+1]-downshift
  }
  text(x2+rightshift,textheight,orderednames,adj=0,cex=0.6)
  
  output<-list(IVnames=IVnames,AUCs=AUCs)
  print(output)
}

#Making the plot
#pdf("Round4_Figure3.pdf",width=7.5,height=7.5)
par(mfrow=c(2,2),mar=c(0.1,5.1,2.1,1.1),mgp=c(4,1,0),las=1, bty="n", xaxt="n")

#########
#Model 1
# start with the full model
m1.AUC.main<-somers2(plogis(predict(m1, newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
plot(0:1,0:1,type="n",xaxt="n",ylim=c(0.55,0.70),xlab="",ylab="Predictive Power (AUC)",main="M1. IDA Project Model")

# Stage 1:
text(0.1,0.70,"Stage 1",adj=0)
points(0,m1.AUC.main,pch=19)
stage1<-treeplot2(model= m1, DV=data$result,IVs=subset(data, select=c(ida.prop, NR.total.rents.ay.lag1, approval_year.count, env.min.combined, env.sector)),
                  IVnames=c("IDA Project", "Resources", "Year", "Env Ministry", "Env Sector"), newdv=data2$result, olddat=data, newdat=data2, x1=0,y1=m1.AUC.main,x2=0.1,downshift=0.004, rightshift=0.02)

# Stage 2A: (without Environmental Sector)
text(0.45,0.70,"Stage 2",adj=0)
points(0.37,0.6314081,pch=19)
stage1<-treeplot2(model= m1 ,DV=data$result,IVs=subset(data, select=c(ida.prop, NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("IDA Project", "Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.37,y1=0.6314081,x2=0.47,downshift=0.004, rightshift=0.02)

# Stage 3A: (without IDA Proportion or Environmental Sector)
text(0.8,0.70,"Stage 3",adj=0)
points(0.74,0.6196040,pch=19)
stage1<-treeplot2(model= m1 ,DV=data$result,IVs=subset(data, select=c(NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.74,y1=0.6196040,x2=.84,downshift=0.004, rightshift=0.02)

# Now add in minimal model:
model.min<-glmer(result ~ (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
AUC.min<-somers2(plogis(predict(model.min,newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
abline(h=AUC.min,lty=3)
text(0,AUC.min+0.005,"Minimal Model: Country Random Effects Only",adj=0,cex=0.8)

########
#Model 2

# start with the full model
m2.AUC.main<-somers2(plogis(predict(m2, newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
plot(0:1,0:1,type="n",xaxt="n",ylim=c(0.55,0.70),xlab="",ylab="Predictive Power (AUC)",main="M2. IDA Portfolio Model")

# Stage 1:
text(0.1,0.70,"Stage 1",adj=0)
points(0,m2.AUC.main,pch=19)
stage1<-treeplot2(model= m2, DV=data$result,IVs=subset(data, select=c(ida.prop.portfolio.ay.proj, NR.total.rents.ay.lag1, approval_year.count, env.min.combined, env.sector)),
                  IVnames=c("IDA Portfolio", "Resources", "Year", "Env Ministry", "Env Sector"), newdv=data2$result, olddat=data, newdat=data2, x1=0,y1=m2.AUC.main,x2=0.1,downshift=0.004, rightshift=0.02)

# Stage 2A: (without Environmental Sector)
text(0.45,0.70,"Stage 2",adj=0)
points(0.37,0.6275468,pch=19)
stage1<-treeplot2(model= m2 ,DV=data$result,IVs=subset(data, select=c(ida.prop.portfolio.ay.proj, NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("IDA Portfolio", "Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.37,y1=0.6275468,x2=0.47,downshift=0.004, rightshift=0.02)

# Stage 3A: (without IDA Proportion or Environmental Sector)
text(0.8,0.70,"Stage 3",adj=0)
points(0.74,0.6196040,pch=19)
stage1<-treeplot2(model= m2 ,DV=data$result,IVs=subset(data, select=c(NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.74,y1=0.6196040,x2=.84,downshift=0.004, rightshift=0.02)

# Now add in the ida.prop.portfolio minimal model:
model.min<-glmer(result ~ (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
AUC.min<-somers2(plogis(predict(model.min,newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
abline(h=AUC.min,lty=3)
text(0,AUC.min+0.005,"Minimal Model: Country Random Effects Only",adj=0,cex=0.8)

########
#Model 3

# start with the full model
m3.AUC.main<-somers2(plogis(predict(m3, newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
plot(0:1,0:1,type="n",xaxt="n",ylim=c(0.55,0.70),xlab="",ylab="Predictive Power (AUC)",main="M3. IDA-Only Category Model")

# Stage 1:
text(0.1,0.70,"Stage 1",adj=0)
points(0,m3.AUC.main,pch=19)
stage1<-treeplot2(model= m3, DV=data$result,IVs=subset(data, select=c(GDP.k.pc.ay.lag1, NR.total.rents.ay.lag1, approval_year.count, env.min.combined, env.sector)),
                  IVnames=c("IDA-Only", "Resources", "Year", "Env Ministry", "Env Sector"), newdv=data2$result, olddat=data, newdat=data2, x1=0,y1=m3.AUC.main,x2=0.1,downshift=0.004, rightshift=0.02)

# Stage 2A: (without Environmental Sector)
text(0.45,0.70,"Stage 2",adj=0)
points(0.37,0.6299230,pch=19)
stage1<-treeplot2(model= m3 ,DV=data$result,IVs=subset(data, select=c(GDP.k.pc.ay.lag1, NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("IDA-Only", "Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.37,y1=0.6299230,x2=0.47,downshift=0.004, rightshift=0.02)

# Stage 3A: (without GDP per Capita or Environmental Sector)
text(0.8,0.70,"Stage 3",adj=0)
points(0.74,0.6196040  ,pch=19)
stage1<-treeplot2(model= m3 ,DV=data$result,IVs=subset(data, select=c(NR.total.rents.ay.lag1, approval_year.count, env.min.combined)),
                  IVnames=c("Resources", "Year", "Env Ministry"), newdv=data2$result, olddat=data, newdat=data2, x1=.74,y1=0.6196040 ,x2=.84,downshift=0.004, rightshift=0.02)

# Now add in the GDP per Capita minimal model:
model.min<-glmer(result ~ (1|country), family=binomial, data=data, control=glmerControl(optimizer="bobyqa"))
AUC.min<-somers2(plogis(predict(model.min,newdata = data2,type = "response", allow.new.levels = TRUE)),data2$result)[1]
abline(h=AUC.min,lty=3)
text(0,AUC.min+0.005,"Minimal Model: Country Random Effects Only",adj=0,cex=0.8)

#dev.off()

########################################################################
###Figure A1
########################################################################
#pdf("Round2_Figure6.pdf",width=6.5,height=4)
par(mfrow=c(1,2),mar=c(2,3,2,1),mgp=c(4,1,0),las=1, bty="n")

###Barplot of achievement by indicator type
data$far <- ordered(data$far, levels=c("R","A","F"))
data$achieve.char <- ifelse(data$achieve==1, "yes","no")
data$achieve.char <- ifelse(is.na(data$achieve.char),"NA",data$achieve.char)
data$achieve.char <- ordered(data$achieve.char, levels=c("yes","no","NA"))
counts <- table(data$achieve.char,data$far)
barplot(counts, main="Achievement", legend=(c("achieved","not achieved","no baseline")), names.arg=c("Function","Action","Form"),ylim=c(0,400),args.legend=list(cex=1,x=3.7,y=400))

###Barplot of persistence by indicator type
data$far <- ordered(data$far, levels=c("R","A","F"))
data$persist.char <- ifelse(data$persist==1, "yes","no")
data$persist.char <- ifelse(is.na(data$persist.char),"NA",data$persist.char)
data$persist.char <- ordered(data$persist.char, levels=c("yes","no","NA"))
counts2 <- table(data$persist.char,data$far)
barplot(counts2, main="Persistence", legend=(c("persistence","no persistence","not sourced")), names.arg=c("Function","Action","Form"),ylim=c(0,400),args.legend=list(cex=1,x=3.7,y=400))

#dev.off()

########################################################################
###Table A1
########################################################################
###Models of Achievement
mod1 <- glmer(achieve.pos ~ form + action + (1|WB_projectid) + (1|country) + (1|level), family=binomial, data=data)
summary(mod1)

data.mod2 <- na.omit(subset(data, select=c("achieve.pos","form", "action", "es.borrower.sat","WB_projectid","country","level")))
mod2 <- glmer(achieve.pos ~ form + action + es.borrower.sat + (1|WB_projectid) + (1|country) + (1|level), family=binomial, data=data.mod2)
summary(mod2)

data.a <- data[!is.na(data$achieve.pos),]
L <- nrow(data.a)
for (i in 1:L){
  sub <- subset(data.a, WB_projectid==data.a$WB_projectid[i])
  data.a$far.proj.var[i] <- ifelse(length(unique(sub$far))>1,1,0)
}
data.fe.a <- subset(data.a, far.proj.var==1)
mod3 <- glm(achieve.pos ~ form + action + WB_projectid, family=binomial(link="logit"), data=data.fe.a)
summary(mod3)

###Models of Persistence
mod4 <- glmer(persist ~ form + action + (1|WB_projectid) + (1|country) + (1|level), family=binomial, data=data)
summary(mod4)

data.mod5 <- na.omit(subset(data, select=c("persist","form", "action", "es.borrower.sat", "gap", "achieve.pos", "WB_projectid","country","level")))
mod5 <- glmer(persist ~ form + action + es.borrower.sat + gap + achieve.pos + (1|WB_projectid) + (1|country) + (1|level), family=binomial, data=data.mod5)
summary(mod5)

data.p <- data[!is.na(data$persist),]
L <- nrow(data.p)
for (i in 1:L){
  sub <- subset(data.p, WB_projectid==data.p$WB_projectid[i])
  data.p$far.proj.var[i] <- ifelse(length(unique(sub$far))>1,1,0)
}
data.fe.p <- subset(data.p, far.proj.var==1)
mod6 <- glm(persist ~ form + action + WB_projectid, family=binomial(link="logit"), data=data.fe.p)
summary(mod6)

########################################################################
###Figure A2
########################################################################
sim.gate <- sim(mod1,n.sims=1000)
sim.gate2 <- sim(mod4,n.sims=1000)

#pdf(file="AchieveMaintain_FirstDifferences.pdf", width=6.5,height=5)
par(mfrow=c(1,2),oma=c(0,0,0,0.1),mar=c(4.3,3.8,3,0))

fixef <- fixef(sim.gate)
ranef.project <- rnorm(1000,mean=0,sd=1.08296) #Taken from model print-out
ranef.country <- rnorm(1000,mean=0,sd=0.49035) #Taken from model print-out
sim.result <- invlogit(ranef.project + ranef.country + fixef[,1] + fixef[,2]*0 + fixef[,3]*0)
sim.form <- invlogit(ranef.project + ranef.country + fixef[,1] + fixef[,2]*1 + fixef[,3]*0)
sim.diff.achieve <- sim.form-sim.result
plot(density(sim.diff.achieve),lwd=3,xlim=c(-0.1,0.35),
     xlab="Change in Probability of Achievement \nby Choosing Form Target", cex.lab=0.8,
     main="Average baseline probability of \nachieving result target = 0.73",cex.main=0.8)
abline(v=0,lty=2,lwd=5,col="red")
text(x=-0.06,y=7,paste("6.5%"))
text(x=0.13,y=7,paste("93.5%"))

fixef2 <- fixef(sim.gate2)
ranef.project2 <- rnorm(1000,mean=0,sd=0.31453) #Taken from model print-out
ranef.country2 <- rnorm(1000,mean=0,sd=0.0000028441) #Taken from model print-out
sim.result2 <- invlogit(ranef.project2 + ranef.country2 + fixef2[,1] + fixef2[,2]*0 + fixef2[,3]*0)
sim.form2 <- invlogit(ranef.project2 + ranef.country2 + fixef2[,1] + fixef2[,2]*1 + fixef2[,3]*0)
sim.diff.maintain <- sim.form2-sim.result2
plot(density(sim.diff.maintain),lwd=3,xlim=c(-0.05,0.7),
     xlab="Change in Probability of Persistence \nby Choosing Form Target",cex.lab=0.8,
     main="Average baseline probability of \nmaintaining result target = 0.59",cex.main=0.8,
     ylab="")
abline(v=0,lty=2,lwd=5,col="red")
text(x=0.6,y=4,paste("100%"))
#dev.off()

########################################################################
###Table B1
########################################################################

sub <- subset(data,far!="A") #Removing "action" targets from the sample
sub$NR.total.rents.ay.lag1.prop <- sub$NR.total.rents.ay.lag1/100

m1.b1 <- glmer(result ~ ida.prop + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
summary(m1.b1)

m2.b1 <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1.prop + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
#Note: NR.total.rents.ay.lag1 rescaled for covergence, but reported on original scale for consistency with rest of paper/appendix
summary(m2.b1)

m3.b1 <- glmer(result ~ op.cat.I + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
summary(m3.b1)

########################################################################
###Figure D1
########################################################################
band.test <- data.frame(seq(from=100, to=500, by=5))
names(band.test) <- "bandwidth"

for (i in 1:nrow(band.test)){
  
  sub <- subset(data, abs(cutoff.diff)<=band.test$bandwidth[i])
  band.test$n[i] <- nrow(sub)
  
  m1.band <- NA
  m1.band <- glmer(result ~ ida.prop + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
  band.test$ida.prop[i]<-fixef(m1.band)[2]
  band.test$ida.prop.se[i]<-se.fixef(m1.band)[2]
  
  m2.band <- NA
  m2.band <- glmer(result ~ ida.prop.portfolio.ay.proj + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
  band.test$ida.prop.portfolio.ay.proj[i]<-fixef(m2.band)[2]
  band.test$ida.prop.portfolio.ay.proj.se[i]<-se.fixef(m2.band)[2]
  
  m3.band <- NA
  m3.band <- glmer(result ~ op.cat.I + NR.total.rents.ay.lag1 + env.min.combined + approval_year.count + env.sector + (1|level) + (1|country), family=binomial, data=sub, control=glmerControl(optimizer="bobyqa"))
  band.test$op.cat.I[i]<-fixef(m3.band)[2]
  band.test$op.cat.I.se[i]<-se.fixef(m3.band)[2]
  
}

jo <- loess(band.test$ida.prop~band.test$bandwidth)
jo.se.high <- loess((band.test$ida.prop + 1.645*band.test$ida.prop.se)~band.test$bandwidth)
jo.se.low <- loess((band.test$ida.prop - 1.645*band.test$ida.prop.se)~band.test$bandwidth)

plot(x=band.test$bandwidth, y=band.test$ida.prop, ylim=c(-5.5,1), xlab="Distance from IDA-Only Eligibility Cutoff (Bandwidth, GDP per capita $)",
     ylab="Coefficient Estimate", main="IDA Proportion at Project Level")
lines(y=predict(jo),x=band.test$bandwidth,col="red", lwd=3)
lines(y=predict(jo.se.high),x=band.test$bandwidth,col="light blue", lwd=3)
lines(y=predict(jo.se.low),x=band.test$bandwidth,col="light blue", lwd=3)
abline(h=0)

ko <- loess(band.test$ida.prop.portfolio.ay.proj~band.test$bandwidth)
ko.se.high <- loess((band.test$ida.prop.portfolio.ay.proj + 1.645*band.test$ida.prop.portfolio.ay.proj.se)~band.test$bandwidth)
ko.se.low <- loess((band.test$ida.prop.portfolio.ay.proj - 1.645*band.test$ida.prop.portfolio.ay.proj.se)~band.test$bandwidth)

plot(x=band.test$bandwidth, y=band.test$ida.prop.portfolio.ay.proj, ylim=c(-6.5,2), xlab="Distance from IDA-Only Eligibility Cutoff (Bandwidth)",
     ylab="Coefficient Estimate", main="IDA Proportion in Country Portfolio")
lines(y=predict(ko),x=band.test$bandwidth,col="red", lwd=3)
lines(y=predict(ko.se.high),x=band.test$bandwidth,col="light blue", lwd=3)
lines(y=predict(ko.se.low),x=band.test$bandwidth,col="light blue", lwd=3)
abline(h=0)

lo <- loess(band.test$op.cat.I~band.test$bandwidth)
lo.se.high <- loess((band.test$op.cat.I + 1.645*band.test$op.cat.I.se)~band.test$bandwidth)
lo.se.low <- loess((band.test$op.cat.I - 1.645*band.test$op.cat.I.se)~band.test$bandwidth)

plot(x=band.test$bandwidth, y=band.test$op.cat.I, ylim=c(-4.5,.5), xlab="Distance from IDA-Only Eligibility Cutoff (Bandwidth)",
     ylab="Coefficient Estimate", main="IDA-Only Categorization")
lines(y=predict(lo),x=band.test$bandwidth,col="red", lwd=3)
lines(y=predict(lo.se.high),x=band.test$bandwidth,col="light blue", lwd=3)
lines(y=predict(lo.se.low),x=band.test$bandwidth,col="light blue", lwd=3)
abline(h=0)
